In [1]:
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from sklearn.model_selection import train_test_split, KFold
from sklearn.metrics import mean_squared_error, r2_score
from sklearn import linear_model
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import median_absolute_error
from matplotlib.pylab import rcParams
import seaborn as sns # more plots
sns.set()
from dateutil.relativedelta import relativedelta # working with dates with style
from scipy.optimize import minimize # for function minimization
import statsmodels.tsa.api as smt
import statsmodels.api as sm
import scipy.stats as scs
from itertools import product # some useful functions
from tqdm import tqdm_notebook
In [2]:
from sqlalchemy import create_engine, text
from datetime import datetime as dt
# Data Corabastos
engine = create_engine('postgres://postgres:postADS4#@192.99.145.205:54328/ADR_DS4A')
precios_df = pd.read_sql_table("adr_corabastos_prices",con=engine,
parse_dates=['cor_date'])
precios_df.columns = ['id', 'Fecha' ,'Categorias','Producto','Presentación','Cantidad','Unidad','Precio Calidad Extra','Precio Calidad Primera','Precio Calidad Corriente','Precio Grandes Superficies']
precios_df['Mes'] = pd.DatetimeIndex(precios_df['Fecha']).month
precios_df['Day'] = pd.DatetimeIndex(precios_df['Fecha']).dayofyear
precios_df['Day_m'] = pd.DatetimeIndex(precios_df['Fecha']).day
precios_df['Period'] = pd.to_datetime(precios_df['Fecha']).dt.to_period("M")
ipc = pd.read_sql_table( "adr_ipc",con=engine)
ipc['Period'] = ipc['ipc_fecha'].dt.to_period("M")
current_ipc = float(ipc[ipc['ipc_fecha']==dt.today().strftime("20%y-%m-01")]['ipc_valor'])
period = str(ipc[ipc['ipc_valor']==current_ipc]['Period'].iloc[0])
df_corriente = precios_df.merge(ipc, on="Period", how='left')
df_corriente['Precio (pesos corrientes)'] = df_corriente['Precio Calidad Primera']/df_corriente['ipc_valor']*current_ipc
In [3]:
###Import Data Clean###
# Price Data
Price_data = pd.read_csv('General/Price_data.csv', index_col=['Year'],
parse_dates=['Year'])
Price_data = Price_data.drop('CU', axis=1)
# Production_data
Production_data=pd.read_csv('General/Production_data.csv', index_col=['Year'],
parse_dates=['Year'])
#Production_Resume_data
Production_Resume_data=pd.read_csv('General/Production_Resume_data.csv', index_col=['Year'],
parse_dates=['Year'])
Top_product=['Leche', 'Leche no procesada','mantequilla y leche','Leche procesada','Leche deslactosada','Leche deslactosada','Leche entera','Leche en polvo','Leche evaporada'
,'Cafe y productos','Cafe','Platanos y otros','Platanos','Yuca','Granos de cacao''Miel''Verduras','Verduras frescas','Verduras de primera','Legumbres',
'Maiz','Aceite de maiz','Aguacate','Arroz']
#Price_data_top=Price_data[Price_data['Item'].str.contains("Cafe")]
Price_data_top=Price_data[Price_data['Item'].isin(Top_product)]
Production_data_top=Production_data[Production_data['Item'].isin(Top_product)]
Production_Resume_data_top=Production_Resume_data[Production_Resume_data['Item'].isin(Top_product)]
In [4]:
g = sns.FacetGrid(Price_data_top,col='Item',col_wrap=2, height=5, despine=True)
g = (g.map(sns.lineplot,data=Price_data_top["Value"],marker=".").set_titles("{col_name} "))
In [5]:
g = sns.FacetGrid(Production_data_top,col='Item',col_wrap=2, height=5, despine=True)
g = (g.map(sns.lineplot,data=Production_data_top["Value"],marker=".").set_titles("{col_name} "))
In [6]:
g = sns.FacetGrid(Production_Resume_data_top,col='Item',col_wrap=2, height=5, despine=True)
g = (g.map(sns.lineplot,data=Production_Resume_data_top["Value"],marker=".").set_titles("{col_name} "))
In [7]:
producto=['PAPA PASTUSA']
precios_df_P = precios_df[precios_df['Producto'].isin(producto)]
df_corriente_P = df_corriente[precios_df['Producto'].isin(producto)]
#df_corriente_P['Fecha']=pd.to_datetime(df_corriente_P['Fecha'])
In [8]:
df_corriente.dtypes
Out[8]:
In [9]:
corriente_Y = df_corriente_P.groupby(df_corriente['Fecha'].dt.to_period('Y')).mean().reset_index()
corriente_M = df_corriente_P.groupby(df_corriente['Fecha'].dt.to_period('M')).mean().reset_index()
corriente_W = df_corriente_P.groupby(df_corriente['Fecha'].dt.to_period('W')).mean().reset_index()
In [10]:
precios_Y = precios_df_P.groupby(precios_df_P['Fecha'].dt.to_period('Y')).mean().reset_index()
precios_M = precios_df_P.groupby(precios_df_P['Fecha'].dt.to_period('M')).mean().reset_index()
precios_W = precios_df_P.groupby(precios_df_P['Fecha'].dt.to_period('W')).mean().reset_index()
In [11]:
precios_Y.head()
Out[11]:
In [12]:
precios_Y['Year'] = precios_Y['Fecha'].dt.year
precios_M['Year'] = precios_M['Fecha'].dt.year
precios_M['Month'] = precios_M['Fecha'].dt.month
precios_M['Theta'] = precios_M['Month']*360/12
precios_W['Month'] = precios_W['Fecha'].dt.month
precios_W['Year'] = precios_W['Fecha'].dt.year
In [13]:
precios_M.columns
Out[13]:
In [14]:
corriente_Y['Year'] = corriente_Y['Fecha'].dt.year
corriente_M['Year'] = corriente_M['Fecha'].dt.year
corriente_M['Month'] = corriente_M['Fecha'].dt.month
corriente_M['Theta'] = corriente_M['Month']*360/12
corriente_W['Month'] = corriente_W['Fecha'].dt.month
corriente_W['Year'] = corriente_W['Fecha'].dt.year
In [15]:
corriente_M.columns
Out[15]:
In [16]:
import plotly.express as px
df = precios_M
fig = px.line_polar(df, r="Precio Calidad Primera", theta="Theta", color="Year", line_close=True,
color_discrete_sequence=px.colors.sequential.Plasma_r,
template="plotly_dark")
fig.show()
In [17]:
import plotly.express as px
df = corriente_M
fig = px.line_polar(df, r="Precio (pesos corrientes)", theta="Theta", color="Year", line_close=True,
color_discrete_sequence = px.colors.sequential.Plasma_r, template="plotly_dark")
fig.show()
In [18]:
### Filtrar un solo producto ###
#Product=['Manzanas']
#Price_data_P=Price_data[Price_data['Item'].isin(Product)].sort_index()
#Production_data_P=Production_data[Production_data['Item'].isin(Product)].sort_index()
#Production_Resume_data_P=Production_Resume_data[Production_Resume_data['Item'].isin(['Tomate'])].sort_index()
#precios_df_P = precios_df[precios_df['Producto'].isin(['CAFE 500GR'])].sort_index()
#corriente_M = df_corriente[df_corriente['Producto'].isin(['CAFE 500GR'])].sort_index()
#Price_data_P=Price_data_P.groupby(Price_data_P.index).mean().rename(columns={0:'mean'})
Moving average¶
In [19]:
corriente_Y['Fecha'] = corriente_Y['Fecha'].astype(str)
corriente_Y['Fecha'] = pd.to_datetime(corriente_Y['Fecha'])
corriente_M['Fecha'] = corriente_M['Fecha'].astype(str)
corriente_M['Fecha'] = pd.to_datetime(corriente_M['Fecha'])
corriente_W['Fecha'] = corriente_W['Fecha'].astype(str).str.split('/',1,expand=True)
corriente_M['Fecha'] = pd.to_datetime(corriente_M['Fecha'])
In [20]:
df_corriente_P.head()
Out[20]:
In [21]:
corriente_Y = corriente_Y.set_index('Fecha')
corriente_M = corriente_M.set_index('Fecha')
corriente_W = corriente_W.set_index('Fecha')
df_corriente_P = df_corriente_P.set_index('Fecha').sort_index()
In [22]:
df_corriente_P.columns
Out[22]:
In [23]:
corriente_Y = corriente_Y .drop([ 'id', 'Cantidad', 'Precio Calidad Extra', 'Precio Calidad Primera',
'Precio Calidad Corriente', 'Precio Grandes Superficies', 'Mes', 'Day',
'Day_m', 'ipc_id', 'ipc_valor', 'Year'], axis=1)
corriente_M = corriente_M.drop(['id', 'Cantidad', 'Precio Calidad Extra', 'Precio Calidad Primera',
'Precio Calidad Corriente', 'Precio Grandes Superficies', 'Mes', 'Day',
'Day_m', 'ipc_id', 'ipc_valor', 'Year',
'Month', 'Theta'], axis=1)
corriente_W = corriente_W.drop(['id', 'Cantidad', 'Precio Calidad Extra', 'Precio Calidad Primera',
'Precio Calidad Corriente', 'Precio Grandes Superficies', 'Mes', 'Day',
'Day_m', 'ipc_id', 'ipc_valor', 'Month','Year'], axis=1)
df_corriente_P = df_corriente_P.drop(['id', 'Categorias', 'Producto', 'Presentación', 'Cantidad', 'Unidad',
'Precio Calidad Extra', 'Precio Calidad Primera',
'Precio Calidad Corriente', 'Precio Grandes Superficies', 'Mes', 'Day',
'Day_m', 'Period', 'ipc_id', 'ipc_fecha', 'ipc_valor'], axis=1)
In [24]:
corriente_Y['Precio (pesos corrientes)']=corriente_Y['Precio (pesos corrientes)'].round(decimals=0).astype(int)
corriente_M['Precio (pesos corrientes)']=corriente_M['Precio (pesos corrientes)'].round(decimals=0).astype(int)
corriente_W['Precio (pesos corrientes)']=corriente_W['Precio (pesos corrientes)'].round(decimals=0).astype(int)
In [25]:
#Price_data_P = Price_data_P.drop(['Item Code'], axis=1)
#precios_df_P = precios_df_P.drop([ 'id', 'Categorias', 'Producto', 'Presentación', 'Cantidad', 'Unidad',
#'Precio Calidad Extra', 'Precio Calidad Corriente', 'Precio Grandes Superficies','Mes','Day','Day_m','Period'], axis=1)
#precios_df_P = precios_df_P.set_index('Fecha')
#precios_df_P = precios_df_P.sort_index(axis = 0)
In [26]:
def plotMovingAverage(series, window, plot_intervals=False, scale=1.96,
plot_anomalies=False):
"""
series - dataframe with timeseries
window - rolling window size
plot_intervals - show confidence intervals
plot_anomalies - show anomalies
"""
rolling_mean = series.rolling(window=window).mean()
plt.figure(figsize=(15,5))
plt.title("Moving average\n window size = {}".format(window))
plt.plot(rolling_mean, "g", label="Rolling mean trend")
# Plot confidence intervals for smoothed values
if plot_intervals:
mae = mean_absolute_error(series[window:], rolling_mean[window:])
deviation = np.std(series[window:] - rolling_mean[window:])
lower_bond = rolling_mean - (mae + scale * deviation)
upper_bond = rolling_mean + (mae + scale * deviation)
plt.plot(upper_bond, "r--", label="Upper Bond / Lower Bond")
plt.plot(lower_bond, "r--")
# Having the intervals, find abnormal values
if plot_anomalies:
anomalies = pd.DataFrame(index=series.index, columns=series.columns)
anomalies[series<lower_bond] = series[series<lower_bond]
anomalies[series>upper_bond] = series[series>upper_bond]
plt.plot(anomalies, "ro", markersize=10)
plt.plot(series[window:], label="Actual values")
plt.legend(loc="upper left")
plt.grid(True)
In [27]:
plotMovingAverage(corriente_Y, 5)
In [28]:
plotMovingAverage(corriente_M, 12)
In [29]:
plotMovingAverage(corriente_W, 5)
In [30]:
plotMovingAverage(df_corriente_P, 7)
In [31]:
plotMovingAverage(corriente_Y, 2, plot_intervals=True,plot_anomalies=True)
In [32]:
plotMovingAverage(corriente_M, 12, plot_intervals=True,plot_anomalies=True)
In [33]:
plotMovingAverage(corriente_W, 7, plot_intervals=True,plot_anomalies=True)
In [34]:
plotMovingAverage(df_corriente_P, 7, plot_intervals=True, plot_anomalies=True)
In [35]:
def plotProcess(x,rho=0):
w = x
for t in range(len(x)):
x[t] = rho * x[t-1] + w[t]
with plt.style.context('bmh'):
plt.figure(figsize=(10, 3))
plt.plot(x)
plt.title("Rho {}\n Dickey-Fuller p-value: {}".format(rho, round(sm.tsa.stattools.adfuller(x)[1], 3)))
In [36]:
for rho in [0, 0.6, 0.9, 1]:
plotProcess(corriente_Y['Precio (pesos corrientes)'].values,rho=rho)
In [37]:
for rho in [0, 0.6, 0.9, 1]:
plotProcess( corriente_M['Precio (pesos corrientes)'].values,rho=rho)
In [38]:
for rho in [0, 0.6, 0.9, 1]:
plotProcess( corriente_W['Precio (pesos corrientes)'].values,rho=rho)
In [39]:
for rho in [0, 0.6, 0.9, 1]:
plotProcess( corriente_W['Precio (pesos corrientes)'].values,rho=rho)
In [40]:
for rho in [0, 0.6, 0.9, 1]:
plotProcess( df_corriente_P['Precio (pesos corrientes)'].values,rho=rho)
Augmented Dickey-Fuller test¶
In [41]:
X = corriente_Y['Precio (pesos corrientes)']
result = sm.tsa.stattools.adfuller(X)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\t%s: %.3f' % (key, value))
In [42]:
X =corriente_M['Precio (pesos corrientes)']
result = sm.tsa.stattools.adfuller(X)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\t%s: %.3f' % (key, value))
In [43]:
X =corriente_W['Precio (pesos corrientes)']
result = sm.tsa.stattools.adfuller(X)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\t%s: %.3f' % (key, value))
In [44]:
X =df_corriente_P['Precio (pesos corrientes)']
result = sm.tsa.stattools.adfuller(X)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\t%s: %.3f' % (key, value))
In [45]:
def tsplot(y, lags=None, figsize=(12, 7), style='bmh'):
"""
Plot time series, its ACF and PACF, calculate Dickey–Fuller test
y - timeseries
lags - how many lags to include in ACF, PACF calculation
"""
if not isinstance(y, pd.Series):
y = pd.Series(y)
with plt.style.context(style):
fig = plt.figure(figsize=figsize)
layout = (2, 2)
ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
acf_ax = plt.subplot2grid(layout, (1, 0))
pacf_ax = plt.subplot2grid(layout, (1, 1))
y.plot(ax=ts_ax)
p_value = sm.tsa.stattools.adfuller(y)[1]
ts_ax.set_title('Time Series Analysis Plots\n Dickey-Fuller: p={0:.5f}'.format(p_value))
smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)
smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)
plt.tight_layout()
In [46]:
tsplot(df_corriente_P['Precio (pesos corrientes)'], lags=6)
In [47]:
tsplot(corriente_M['Precio (pesos corrientes)'], lags=100)
In [48]:
tsplot(corriente_W['Precio (pesos corrientes)'], lags=200)
In [49]:
tsplot(df_corriente_P['Precio (pesos corrientes)'], lags=200)
In [50]:
corriente_Y_diff = corriente_Y['Precio (pesos corrientes)'] - corriente_Y['Precio (pesos corrientes)'].shift(1)
tsplot(corriente_Y_diff[1:], lags=6)
In [51]:
corriente_M_diff = corriente_M['Precio (pesos corrientes)'] - corriente_M['Precio (pesos corrientes)'].shift(40)
tsplot(corriente_M_diff [40:], lags=60)
In [52]:
corriente_W_diff = corriente_W['Precio (pesos corrientes)'] - corriente_W['Precio (pesos corrientes)'].shift(1)
tsplot(corriente_W_diff[1:], lags=60)
In [ ]:
In [53]:
df_corriente_P_diff = df_corriente_P['Precio (pesos corrientes)'] - df_corriente_P['Precio (pesos corrientes)'].shift(180)
tsplot(df_corriente_P_diff[180:], lags=60)
In [54]:
ps = range(1, 4)
d=1
qs = range(1, 4)
Ps = range(0, 3)
D=1
Qs = range(0, 3)
s = 8
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[54]:
In [55]:
def optimizeSARIMA(parameters_list, d, D, s):
"""
Return dataframe with parameters and corresponding AIC
parameters_list - list with (p, q, P, Q) tuples
d - integration order in ARIMA model
D - seasonal integration order
s - length of season
"""
results = []
best_aic = float("inf")
for param in tqdm_notebook(parameters_list):
# we need try-except because on some combinations model fails to converge
try:
model=sm.tsa.statespace.SARIMAX(corriente_M['Precio (pesos corrientes)'], order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], s)).fit(disp=-1)
except:
continue
aic = model.aic
# saving best model, AIC and parameters
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
# sorting in ascending order, the lower AIC is - the better
result_table = result_table.sort_values(by='aic', ascending=True).reset_index(drop=True)
return result_table
In [56]:
%%time
result_table_M = optimizeSARIMA(parameters_list, d, D, s)
In [57]:
result_table_M.head(40)
Out[57]:
In [58]:
p, q, P, Q = result_table_M.parameters[0]
best_model=sm.tsa.statespace.SARIMAX(corriente_M['Precio (pesos corrientes)'], order=(p, d, q),
seasonal_order=(P, D, Q, s)).fit(disp=-1)
print(best_model.summary())
In [59]:
def optimizeSARIMA(parameters_list, d, D, s):
"""
Return dataframe with parameters and corresponding AIC
parameters_list - list with (p, q, P, Q) tuples
d - integration order in ARIMA model
D - seasonal integration order
s - length of season
"""
results = []
best_aic = float("inf")
for param in tqdm_notebook(parameters_list):
# we need try-except because on some combinations model fails to converge
try:
model=sm.tsa.statespace.SARIMAX(corriente_W['Precio (pesos corrientes)'], order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], s)).fit(disp=-1)
except:
continue
aic = model.aic
# saving best model, AIC and parameters
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
# sorting in ascending order, the lower AIC is - the better
result_table = result_table.sort_values(by='aic', ascending=True).reset_index(drop=True)
return result_table
In [ ]:
%%time
result_table_W = optimizeSARIMA(parameters_list, d, D, s)
In [ ]:
result_table_W.head(20)
In [ ]:
p, q, P, Q = result_table_W.parameters[0]
best_model=sm.tsa.statespace.SARIMAX(corriente_W['Precio (pesos corrientes)'], order=(p, d, q),
seasonal_order=(P, D, Q, s)).fit(disp=-1)
print(best_model.summary())
In [ ]:
print(mle_retvals)
In [ ]:
tsplot(best_model.resid[+1:], lags=60)
In [ ]: